
* generate variable medians for NY, used to get predicted probabilities and marginal effects later
 
sum f_c9010t_perpov_dw_z if (ra_site==5 & sample_nox==1), detail
scalar f_c9010t_perpov_dw_z_med_nox = r(p50)

sum f_c9010t_pminorty_dw_z if (ra_site==5 & sample_nox==1), detail
scalar f_c9010t_pminorty_dw_z_med_nox = r(p50)

sum f_c9010t_perpov_dw_z if (ra_site==5 & sample_x==1), detail
scalar f_c9010t_perpov_dw_z_med_x = r(p50)

sum f_c9010t_pminorty_dw_z if (ra_site==5 & sample_x==1), detail
scalar f_c9010t_pminorty_dw_z_med_x = r(p50)

foreach x in $x_covariates {
		
	summarize `x' if (ra_site==5 & sample_x==1), detail
	scalar `x'_med_x = r(p50)
	}
	
sum x_f_release1 if (ra_site==5 & sample_nox==1), detail
scalar x_f_release1_med_nox = r(p50)

sum x_f_release1 if (ra_site==5 & sample_x==1), detail
scalar x_f_release1_med_x = r(p50)

scalar x_rad_ad_41_45_med_x = 1
scalar x_f_hh_size3_med_x = 1

scalar x_rad_ad_ethrace_black_nh_med_x = 0 
scalar x_rad_ad_ethrace_hisp_med_x = 1 

** x_rad_ad_ethrace_black_nh_med = 0 
** x_rad_ad_ethrace_hisp_med = 1 
** tmp_female_med = 1
** x_rad_ad_le_35_med = 0 (9.7% of adults)
** x_rad_ad_36_40_med = 0 (20.7% of adults)
** x_rad_ad_41_45_med = 0 (force set to 1; 23.9% of adults)
** x_rad_ad_46_50_med = 0 (18.1% of adults)
** x_f_ad_nevmarr_med = 1
** x_f_ad_parentu18_med = 0
** x_f_ad_working_med = 0
** x_f_ad_edinsch_med = 0
** x_f_ad_edgradhs_med = 0
** x_f_ad_edged_med = 0
** x_f_hh_afdc_med = 1
** rad_svy_bl_totincm_2009d_med = 11950
** x_f_hh_car_med = 0
** x_f_hh_disabl_med = 0
** x_f_hh_noteens_med = 1
** x_f_hh_size2_med = 0 (22% of adults)
** x_f_hh_size3_med = 0 (force set to 1; 29% of adults)
** x_f_hh_size4_med = 0 (25% of adults)
** x_f_hh_victim_med = 1
** x_f_hood_unsafenit_med = 1
** x_f_hood_verydissat_med = 1
** x_f_hood_5y_med = 1
** x_f_hous_mov3tm_med = 0
** x_f_hood_nofamily_med = 1
** x_f_hood_nofriend_med = 0
** x_f_hood_chat_med = 1
** x_f_hood_nbrkid_med = 0
** x_f_hous_fndapt_med = 0
** x_f_hous_sec8bef_med = 0
** x_f_hous_movdrgs_med = 1
** x_f_hous_movschl_med = 1
** cov_hous_movapt_med = 0
** cov_hous_movjob_med = 0
** x_f_release1_med = 0

* expressions to be used later for predicted probabilities and marginal effects

global x_fix_short "x_f_site_balt=0 x_f_site_bos=0 x_f_site_chi=0 x_f_site_la=0"
global x_fix_long  "x_rad_ad_ethrace_black_nh=0 x_rad_ad_ethrace_hisp=1 tmp_female=1 x_rad_ad_le_35=0 x_rad_ad_36_40=0 x_rad_ad_41_45=1 x_rad_ad_46_50=0 x_f_ad_nevmarr=1 x_f_ad_parentu18=0 x_f_ad_working=0 x_f_ad_edinsch=0 x_f_ad_edgradhs=0 x_f_ad_edged=0 x_f_hh_afdc=1 rad_svy_bl_totincm_2009d=11950 x_f_hh_car=0 x_f_hh_disabl=0 x_f_hh_noteens=1 x_f_hh_size2=0 x_f_hh_size3=1 x_f_hh_size4=0 x_f_hh_victim=1 x_f_hood_unsafenit=1 x_f_hood_verydissat=1 x_f_hood_5y=1 x_f_hous_mov3tm=0 x_f_hood_nofamily=1 x_f_hood_nofriend=0 x_f_hood_chat=1 x_f_hood_nbrkid=0 x_f_hous_fndapt=0 x_f_hous_sec8bef=0 x_f_hous_movdrgs=1 x_f_hous_movschl=1 cov_hous_movapt=0 cov_hous_movjob=0 x_f_release1=0 x_f_site_balt=0 x_f_site_bos=0 x_f_site_chi=0 x_f_site_la=0"

global xgamma_short "(0) "
global xgamma_morex "x_f_release1_med_x*[xb1]x_f_release1+tmp_female_med_x*[xb1]tmp_female+x_rad_ad_le_35_med_x*[xb1]x_rad_ad_le_35+x_rad_ad_36_40_med_x*[xb1]x_rad_ad_36_40+x_rad_ad_41_45_med_x*[xb1]x_rad_ad_41_45+x_rad_ad_46_50_med_x*[xb1]x_rad_ad_46_50+x_rad_ad_ethrace_black_nh_med_x*[xb1]x_rad_ad_ethrace_black_nh+x_rad_ad_ethrace_hisp_med_x*[xb1]x_rad_ad_ethrace_hisp+x_f_ad_nevmarr_med_x*[xb1]x_f_ad_nevmarr+x_f_ad_parentu18_med_x*[xb1]x_f_ad_parentu18+x_f_ad_working_med_x*[xb1]x_f_ad_working+x_f_ad_edinsch_med_x*[xb1]x_f_ad_edinsch+x_f_ad_edgradhs_med_x*[xb1]x_f_ad_edgradhs+x_f_ad_edged_med_x*[xb1]x_f_ad_edged+x_f_hh_afdc_med_x*[xb1]x_f_hh_afdc+rad_svy_bl_totincm_2009d_med_x*[xb1]rad_svy_bl_totincm_2009d+x_f_hh_car_med_x*[xb1]x_f_hh_car+x_f_hh_disabl_med_x*[xb1]x_f_hh_disabl+x_f_hh_noteens_med_x*[xb1]x_f_hh_noteens+x_f_hh_size2_med_x*[xb1]x_f_hh_size2+x_f_hh_size3_med_x*[xb1]x_f_hh_size3+x_f_hh_size4_med_x*[xb1]x_f_hh_size4+x_f_hh_victim_med_x*[xb1]x_f_hh_victim+x_f_hood_unsafenit_med_x*[xb1]x_f_hood_unsafenit+x_f_hood_verydissat_med_x*[xb1]x_f_hood_verydissat+x_f_hood_5y_med_x*[xb1]x_f_hood_5y+x_f_hous_mov3tm_med_x*[xb1]x_f_hous_mov3tm+x_f_hood_nofamily_med_x*[xb1]x_f_hood_nofamily+x_f_hood_nofriend_med_x*[xb1]x_f_hood_nofriend+x_f_hood_chat_med_x*[xb1]x_f_hood_chat+x_f_hood_nbrkid_med_x*[xb1]x_f_hood_nbrkid+x_f_hous_fndapt_med_x*[xb1]x_f_hous_fndapt+x_f_hous_sec8bef_med_x*[xb1]x_f_hous_sec8bef+x_f_hous_movdrgs_med_x*[xb1]x_f_hous_movdrgs+x_f_hous_movschl_med_x*[xb1]x_f_hous_movschl+cov_hous_movapt_med_x*[xb1]cov_hous_movapt+cov_hous_movjob_med_x*[xb1]cov_hous_movjob"

** ordered probit probabilities and marginal effects

** FIGURE 1 

** W = hood poverty, no covariates

** CT model 

ml model lf myoprobit2_lf (xb1:happy_scale012 = $site_covs f_c9010t_perpov_dw_z, nocons) (cut1:) (lndiff:) /*
*/(xb2:f_c9010t_perpov_dw_z = interaction*) (lnsigma:) (rho:) [pw=$wt] , /*
*/diparm(lndiff cut1, f(exp(@1)+@2) d(exp(@1) -exp(@1)+1) label(cut2)) /*
*/diparm(lnsigma, exp label(sigma)) 
ml init xb1:x_f_site_bal=0.1825 x_f_site_bos=0.0528 x_f_site_chi=0.1811 x_f_site_la=0.0674 f_c9010t_perpov_dw_z=-0.1610
ml maximize, difficult
est store pov_t_nox

forvalues y=0/2{
	forvalues i = 1/11{
		scalar poverty_`y'_`i' = `i' - 7
		if `y'==0{
		nlcom normal([cut1]_cons-poverty_`y'_`i'*[xb1]f_c9010t_perpov_dw_z - (${xgamma_short}))
		matrix b = r(b)
		matrix V = r(V)
		scalar py`y'_nox_`i'_med_ny = b[1,1]
		scalar v`y'_nox_`i' = sqrt(V[1,1])
		scalar lci`y'_nox_`i'_med_ny = py`y'_nox_`i'_med_ny-invnormal(0.975)*v`y'_nox_`i'
		scalar uci`y'_nox_`i'_med_ny = py`y'_nox_`i'_med_ny+invnormal(0.975)*v`y'_nox_`i'
		}
		else if `y'==1{
		nlcom normal(exp([lndiff]_cons)+[cut1]_cons-poverty_`y'_`i'*[xb1]f_c9010t_perpov_dw_z - (${xgamma_short}))- /*				
				*/normal([cut1]_cons-poverty_`y'_`i'*[xb1]f_c9010t_perpov_dw_z - (${xgamma_short}))
		matrix b = r(b)
		matrix V = r(V)
		scalar py`y'_nox_`i'_med_ny = b[1,1]
		scalar v`y'_nox_`i' = sqrt(V[1,1])
		scalar lci`y'_nox_`i'_med_ny = py`y'_nox_`i'_med_ny-invnormal(0.975)*v`y'_nox_`i'
		scalar uci`y'_nox_`i'_med_ny = py`y'_nox_`i'_med_ny+invnormal(0.975)*v`y'_nox_`i'
		}
		else if `y'==2{
		nlcom 1-normal(exp([lndiff]_cons)+[cut1]_cons-poverty_`y'_`i'*[xb1]f_c9010t_perpov_dw_z - (${xgamma_short})) 
				
		matrix b = r(b)
		matrix V = r(V)
		scalar py`y'_nox_`i'_med_ny = b[1,1]
		scalar v`y'_nox_`i' = sqrt(V[1,1])
		scalar lci`y'_nox_`i'_med_ny = py`y'_nox_`i'_med_ny-invnormal(0.975)*v`y'_nox_`i'
		scalar uci`y'_nox_`i'_med_ny = py`y'_nox_`i'_med_ny+invnormal(0.975)*v`y'_nox_`i'
		}
	}

	matrix py`y'_nox_pov_med_ny = (scalar(poverty_`y'_1),scalar(py`y'_nox_1_med_ny),scalar(lci`y'_nox_1_med_ny),scalar(uci`y'_nox_1_med_ny)\/*
				*/scalar(poverty_`y'_2),scalar(py`y'_nox_2_med_ny),scalar(lci`y'_nox_2_med_ny),scalar(uci`y'_nox_2_med_ny)\/*
				*/scalar(poverty_`y'_3),scalar(py`y'_nox_3_med_ny),scalar(lci`y'_nox_3_med_ny),scalar(uci`y'_nox_3_med_ny)\/*
				*/scalar(poverty_`y'_4),scalar(py`y'_nox_4_med_ny),scalar(lci`y'_nox_4_med_ny),scalar(uci`y'_nox_4_med_ny)\/*
				*/scalar(poverty_`y'_5),scalar(py`y'_nox_5_med_ny),scalar(lci`y'_nox_5_med_ny),scalar(uci`y'_nox_5_med_ny)\/*
				*/scalar(poverty_`y'_6),scalar(py`y'_nox_6_med_ny),scalar(lci`y'_nox_6_med_ny),scalar(uci`y'_nox_6_med_ny)\/*
				*/scalar(poverty_`y'_7),scalar(py`y'_nox_7_med_ny),scalar(lci`y'_nox_7_med_ny),scalar(uci`y'_nox_7_med_ny)\/*
				*/scalar(poverty_`y'_8),scalar(py`y'_nox_8_med_ny),scalar(lci`y'_nox_8_med_ny),scalar(uci`y'_nox_8_med_ny)\/*
				*/scalar(poverty_`y'_9),scalar(py`y'_nox_9_med_ny),scalar(lci`y'_nox_9_med_ny),scalar(uci`y'_nox_9_med_ny)\/*
				*/scalar(poverty_`y'_10),scalar(py`y'_nox_10_med_ny),scalar(lci`y'_nox_10_med_ny),scalar(uci`y'_nox_10_med_ny)\/*
				*/scalar(poverty_`y'_11),scalar(py`y'_nox_11_med_ny),scalar(lci`y'_nox_11_med_ny),scalar(uci`y'_nox_11_med_ny))

	matrix colnames py`y'_nox_pov_med_ny = pov`y'_nox_med_ny py`y'_nox_med_ny lci`y'_nox_med_ny uci`y'_nox_med_ny
	svmat py`y'_nox_pov_med_ny, names(col)
	
	matrix list py`y'_nox_pov_med_ny
}

** ordered probit model 

ml model lf myoprobit_lf (xb1:happy_scale012 = $site_covs f_c9010t_perpov_dw_z, nocons) (cut1:) (cut2:) [pw=$wt]
ml maximize
est store pov_op_nox

forvalues y = 0/2 {
	if `y'==0{
	margins, at(${x_fix_short} f_c9010t_perpov_dw_z=(-6(1)4)) expression(normal(predict(equation(cut1))-predict(equation(xb1))))	
	}
	else if `y'==1{
	margins, at(${x_fix_short} f_c9010t_perpov_dw_z=(-6(1)4)) expression((normal(predict(equation(cut2))-predict(equation(xb1))))-(normal(predict(equation(cut1))-predict(equation(xb1)))))
	}
	else if `y'==2{
	margins, at(${x_fix_short} f_c9010t_perpov_dw_z=(-6(1)4)) expression(1-(normal(predict(equation(cut2))-predict(equation(xb1)))))
	}
	
	marginsplot, ylabel(-0.4(0.4)1.2, grid) title(" ") xtitle(" ") ytitle(" ") xlabel(-6(2)4) recast(line) recastci(rarea) plotopts(lcolor(black) lwidth(medthick)) ciopts(color(black%50)) ///
	addplot(rarea lci`y'_nox_med_ny uci`y'_nox_med_ny pov`y'_nox_med_ny, ylabel(-0.4(0.4)1.2, grid) color(black%20)||line py`y'_nox_med_ny pov`y'_nox_med_ny, ylabel(-0.4(0.4)1.2, grid) ///
	lcolor(gs5) lpattern(dash) lwidth(medthick)) legend(off)
	graph export ${pathres}prob_y`y'_nox_pov_med_ny_${date}.pdf, replace
	}
	
** W = hood minority, no covariates

** CT model 

ml model lf myoprobit2_lf (xb1:happy_scale012 = $site_covs f_c9010t_pminorty_dw_z, nocons) (cut1:) (lndiff:) /*
*/(xb2:f_c9010t_pminorty_dw_z = interaction*) (lnsigma:) (rho:) [pw=$wt] , /*
*/diparm(lndiff cut1, f(exp(@1)+@2) d(exp(@1) -exp(@1)+1) label(cut2)) /*
*/diparm(lnsigma, exp label(sigma)) 
ml maximize
est store min_t_nox

forvalues y = 0/2{
	forvalues i = 1/11{
		scalar minority_`y'_`i' = `i' - 7
		if `y'==0{
		nlcom normal([cut1]_cons-minority_`y'_`i'*[xb1]f_c9010t_pminorty_dw_z - (${xgamma_short}) )
		}
		else if `y'==1{
		nlcom normal(exp([lndiff]_cons)+[cut1]_cons-minority_`y'_`i'*[xb1]f_c9010t_pminorty_dw_z - (${xgamma_short}) ) - /*				
				*/normal([cut1]_cons-minority_`y'_`i'*[xb1]f_c9010t_pminorty_dw_z - (${xgamma_short}) )
		}
		else if `y'==2{
		nlcom 1-normal(exp([lndiff]_cons)+[cut1]_cons-minority_`y'_`i'*[xb1]f_c9010t_pminorty_dw_z - (${xgamma_short}) ) 
		}
		
		matrix b = r(b)
		matrix V = r(V)
		scalar py`y'_nox_`i'_med_ny = b[1,1]
		scalar v`y'_nox_`i' = sqrt(V[1,1])
		scalar lci`y'_nox_`i'_med_ny = py`y'_nox_`i'_med_ny-invnormal(0.975)*v`y'_nox_`i'
		scalar uci`y'_nox_`i'_med_ny = py`y'_nox_`i'_med_ny+invnormal(0.975)*v`y'_nox_`i'
	}

	matrix py`y'_nox_min_med_ny = (scalar(minority_`y'_1),scalar(py`y'_nox_1_med_ny),scalar(lci`y'_nox_1_med_ny),scalar(uci`y'_nox_1_med_ny)\/*
				*/scalar(minority_`y'_2),scalar(py`y'_nox_2_med_ny),scalar(lci`y'_nox_2_med_ny),scalar(uci`y'_nox_2_med_ny)\/*
				*/scalar(minority_`y'_3),scalar(py`y'_nox_3_med_ny),scalar(lci`y'_nox_3_med_ny),scalar(uci`y'_nox_3_med_ny)\/*
				*/scalar(minority_`y'_4),scalar(py`y'_nox_4_med_ny),scalar(lci`y'_nox_4_med_ny),scalar(uci`y'_nox_4_med_ny)\/*
				*/scalar(minority_`y'_5),scalar(py`y'_nox_5_med_ny),scalar(lci`y'_nox_5_med_ny),scalar(uci`y'_nox_5_med_ny)\/*
				*/scalar(minority_`y'_6),scalar(py`y'_nox_6_med_ny),scalar(lci`y'_nox_6_med_ny),scalar(uci`y'_nox_6_med_ny)\/*
				*/scalar(minority_`y'_7),scalar(py`y'_nox_7_med_ny),scalar(lci`y'_nox_7_med_ny),scalar(uci`y'_nox_7_med_ny)\/*
				*/scalar(minority_`y'_8),scalar(py`y'_nox_8_med_ny),scalar(lci`y'_nox_8_med_ny),scalar(uci`y'_nox_8_med_ny)\/*
				*/scalar(minority_`y'_9),scalar(py`y'_nox_9_med_ny),scalar(lci`y'_nox_9_med_ny),scalar(uci`y'_nox_9_med_ny)\/*
				*/scalar(minority_`y'_10),scalar(py`y'_nox_10_med_ny),scalar(lci`y'_nox_10_med_ny),scalar(uci`y'_nox_10_med_ny)\/*
				*/scalar(minority_`y'_11),scalar(py`y'_nox_11_med_ny),scalar(lci`y'_nox_11_med_ny),scalar(uci`y'_nox_11_med_ny))

	matrix colnames py`y'_nox_min_med_ny = min`y'_nox_med_ny minpy`y'_nox_med_ny minlci`y'_nox_med_ny minuci`y'_nox_med_ny
	svmat py`y'_nox_min_med_ny, names(col)
	
	matrix list py`y'_nox_min_med_ny
}

** ordered probit model 

ml model lf myoprobit_lf (xb1:happy_scale012 = $site_covs f_c9010t_pminorty_dw_z, nocons) (cut1:) (cut2:) [pw=$wt]
ml maximize
est store min_op_nox

forvalues y = 0/2 {
	if `y'==0{
	margins, at(${x_fix_short} f_c9010t_pminorty_dw_z=(-6(1)4)) expression(normal(predict(equation(cut1))-predict(equation(xb1))))
	}
	else if `y'==1{
	margins, at(${x_fix_short} f_c9010t_pminorty_dw_z=(-6(1)4)) expression((normal(predict(equation(cut2))-predict(equation(xb1))))-(normal(predict(equation(cut1))-predict(equation(xb1)))))
	}
	else if `y'==2{
	margins, at(${x_fix_short} f_c9010t_pminorty_dw_z=(-6(1)4)) expression(1-(normal(predict(equation(cut2))-predict(equation(xb1)))))
	}
	
	marginsplot, ylabel(-0.4(0.4)1.2, grid) title(" ") xtitle(" ") ytitle(" ") xlabel(-6(2)4) recast(line) recastci(rarea) plotopts(lcolor(black) lwidth(medthick)) ciopts(color(black%50)) ///
	addplot(rarea minlci`y'_nox_med_ny minuci`y'_nox_med_ny min`y'_nox_med_ny, ylabel(-0.4(0.4)1.2, grid) color(black%20)||line minpy`y'_nox_med_ny min`y'_nox_med_ny, ///
	ylabel(-0.4(0.4)1.2, grid) lcolor(gs5) lpattern(dash) lwidth(medthick)) legend(off)
	graph export ${pathres}prob_y`y'_nox_min_med_ny_${date}.pdf, replace
	}

** FIGURE 2 
	
** W = (hood poverty, hood minority), no covariates

** CT model 

ml model lf myoprobit3_lf (xb1:happy_scale012 = $site_covs f_c9010t_perpov_dw_z f_c9010t_pminorty_dw_z, nocons) (cut1:) (lndiff:) /*
						*/(xb2:f_c9010t_perpov_dw_z = interaction*) (lnvar_v1:) (rho_uv1:) /*
						*/(xb3:f_c9010t_pminorty_dw_z = interaction*) (lnvar_v2:) (rho_uv2:) (cov_v1v2:) [pw=$wt], /*
						*/diparm(lndiff cut1, f(exp(@1)+@2) d(exp(@1) -exp(@1)+1) label(cut2)) /*
						*/diparm(lnvar_v1, exp label(var_v1)) /*
						*/diparm(lnvar_v2, exp label(var_v2)) 
ml init xb1:x_f_site_bal=0.3039 x_f_site_bos=0.3386 x_f_site_chi=0.1527 x_f_site_la=0.0639 f_c9010t_perpov_dw_z=-0.2998 f_c9010t_pminorty_dw_z=0.3151 cut1:_cons=-0.8165
ml maximize, difficult
est store povmin_t_nox

forvalues y = 0/2{
	forvalues i = 1/11{
		scalar poverty_`y'_`i' = `i' - 7
		if `y'==0{
		nlcom normal([cut1]_cons-poverty_`y'_`i'*[xb1]f_c9010t_perpov_dw_z-f_c9010t_pminorty_dw_z_med_nox*[xb1]f_c9010t_pminorty_dw_z-(${xgamma_short}))
		matrix b = r(b)
		matrix V = r(V)
		scalar povmin1_nox_py`y'_`i'_med_ny = b[1,1]
		scalar povmin1_nox_v`y'_`i' = sqrt(V[1,1])
		scalar povmin1_nox_lci`y'_`i'_med_ny = povmin1_nox_py`y'_`i'_med_ny-invnormal(0.975)*povmin1_nox_v`y'_`i'
		scalar povmin1_nox_uci`y'_`i'_med_ny = povmin1_nox_py`y'_`i'_med_ny+invnormal(0.975)*povmin1_nox_v`y'_`i'
		}
		else if `y'==1{
		nlcom normal(exp([lndiff]_cons)+[cut1]_cons-poverty_`y'_`i'*[xb1]f_c9010t_perpov_dw_z-f_c9010t_pminorty_dw_z_med_nox*[xb1]f_c9010t_pminorty_dw_z-(${xgamma_short}))- /*
				*/normal([cut1]_cons-poverty_`y'_`i'*[xb1]f_c9010t_perpov_dw_z-f_c9010t_pminorty_dw_z_med_nox*[xb1]f_c9010t_pminorty_dw_z-(${xgamma_short}))
		matrix b = r(b)
		matrix V = r(V)
		scalar povmin1_nox_py`y'_`i'_med_ny = b[1,1]
		scalar povmin1_nox_v`y'_`i' = sqrt(V[1,1])
		scalar povmin1_nox_lci`y'_`i'_med_ny = povmin1_nox_py`y'_`i'_med_ny-invnormal(0.975)*povmin1_nox_v`y'_`i'
		scalar povmin1_nox_uci`y'_`i'_med_ny = povmin1_nox_py`y'_`i'_med_ny+invnormal(0.975)*povmin1_nox_v`y'_`i'
		}
		else if `y'==2{
		nlcom 1-normal(exp([lndiff]_cons)+[cut1]_cons-poverty_`y'_`i'*[xb1]f_c9010t_perpov_dw_z-f_c9010t_pminorty_dw_z_med_nox*[xb1]f_c9010t_pminorty_dw_z-(${xgamma_short}))
		matrix b = r(b)
		matrix V = r(V)
		scalar povmin1_nox_py`y'_`i'_med_ny = b[1,1]
		scalar povmin1_nox_v`y'_`i' = sqrt(V[1,1])
		scalar povmin1_nox_lci`y'_`i'_med_ny = povmin1_nox_py`y'_`i'_med_ny-invnormal(0.975)*povmin1_nox_v`y'_`i'
		scalar povmin1_nox_uci`y'_`i'_med_ny = povmin1_nox_py`y'_`i'_med_ny+invnormal(0.975)*povmin1_nox_v`y'_`i'
		}
	}

	matrix py`y'_nox_povmin1_med_ny = (scalar(poverty_`y'_1),scalar(povmin1_nox_py`y'_1_med_ny),scalar(povmin1_nox_lci`y'_1_med_ny),scalar(povmin1_nox_uci`y'_1_med_ny)\/*
				*/scalar(poverty_`y'_2),scalar(povmin1_nox_py`y'_2_med_ny),scalar(povmin1_nox_lci`y'_2_med_ny),scalar(povmin1_nox_uci`y'_2_med_ny)\/*
				*/scalar(poverty_`y'_3),scalar(povmin1_nox_py`y'_3_med_ny),scalar(povmin1_nox_lci`y'_3_med_ny),scalar(povmin1_nox_uci`y'_3_med_ny)\/*
				*/scalar(poverty_`y'_4),scalar(povmin1_nox_py`y'_4_med_ny),scalar(povmin1_nox_lci`y'_4_med_ny),scalar(povmin1_nox_uci`y'_4_med_ny)\/*
				*/scalar(poverty_`y'_5),scalar(povmin1_nox_py`y'_5_med_ny),scalar(povmin1_nox_lci`y'_5_med_ny),scalar(povmin1_nox_uci`y'_5_med_ny)\/*
				*/scalar(poverty_`y'_6),scalar(povmin1_nox_py`y'_6_med_ny),scalar(povmin1_nox_lci`y'_6_med_ny),scalar(povmin1_nox_uci`y'_6_med_ny)\/*
				*/scalar(poverty_`y'_7),scalar(povmin1_nox_py`y'_7_med_ny),scalar(povmin1_nox_lci`y'_7_med_ny),scalar(povmin1_nox_uci`y'_7_med_ny)\/*
				*/scalar(poverty_`y'_8),scalar(povmin1_nox_py`y'_8_med_ny),scalar(povmin1_nox_lci`y'_8_med_ny),scalar(povmin1_nox_uci`y'_8_med_ny)\/*
				*/scalar(poverty_`y'_9),scalar(povmin1_nox_py`y'_9_med_ny),scalar(povmin1_nox_lci`y'_9_med_ny),scalar(povmin1_nox_uci`y'_9_med_ny)\/*
				*/scalar(poverty_`y'_10),scalar(povmin1_nox_py`y'_10_med_ny),scalar(povmin1_nox_lci`y'_10_med_ny),scalar(povmin1_nox_uci`y'_10_med_ny)\/*
				*/scalar(poverty_`y'_11),scalar(povmin1_nox_py`y'_11_med_ny),scalar(povmin1_nox_lci`y'_11_med_ny),scalar(povmin1_nox_uci`y'_11_med_ny))

	matrix colnames py`y'_nox_povmin1_med_ny = povmin1_nox_pov`y'_med_ny povmin1_nox_py`y'_med_ny povmin1_nox_lci`y'_med_ny povmin1_nox_uci`y'_med_ny
	svmat py`y'_nox_povmin1_med_ny, names(col)

	matrix list py`y'_nox_povmin1_med_ny
}

forvalues y = 0/2{
	forvalues i = 1/11{
		scalar minority_`y'_`i' = `i' - 7
		if `y'==0{
		nlcom normal([cut1]_cons-f_c9010t_perpov_dw_z_med_nox*[xb1]f_c9010t_perpov_dw_z-minority_`y'_`i'*[xb1]f_c9010t_pminorty_dw_z-(${xgamma_short}))
		matrix b = r(b)
		matrix V = r(V)
		scalar povmin2_nox_py`y'_`i'_med_ny = b[1,1]
		scalar povmin2_nox_v`y'_`i' = sqrt(V[1,1])
		scalar povmin2_nox_lci`y'_`i'_med_ny = povmin2_nox_py`y'_`i'_med_ny-invnormal(0.975)*povmin2_nox_v`y'_`i'
		scalar povmin2_nox_uci`y'_`i'_med_ny = povmin2_nox_py`y'_`i'_med_ny+invnormal(0.975)*povmin2_nox_v`y'_`i'
		}
		else if `y'==1{
		nlcom normal(exp([lndiff]_cons)+[cut1]_cons-f_c9010t_perpov_dw_z_med_nox*[xb1]f_c9010t_perpov_dw_z-minority_`y'_`i'*[xb1]f_c9010t_pminorty_dw_z-(${xgamma_short}))- /*
				*/normal([cut1]_cons-f_c9010t_perpov_dw_z_med_nox*[xb1]f_c9010t_perpov_dw_z-minority_`y'_`i'*[xb1]f_c9010t_pminorty_dw_z-(${xgamma_short}))
		matrix b = r(b)
		matrix V = r(V)
		scalar povmin2_nox_py`y'_`i'_med_ny = b[1,1]
		scalar povmin2_nox_v`y'_`i' = sqrt(V[1,1])
		scalar povmin2_nox_lci`y'_`i'_med_ny = povmin2_nox_py`y'_`i'_med_ny-invnormal(0.975)*povmin2_nox_v`y'_`i'
		scalar povmin2_nox_uci`y'_`i'_med_ny = povmin2_nox_py`y'_`i'_med_ny+invnormal(0.975)*povmin2_nox_v`y'_`i'
		}
		else if `y'==2{
		nlcom 1-normal(exp([lndiff]_cons)+[cut1]_cons-f_c9010t_perpov_dw_z_med_nox*[xb1]f_c9010t_perpov_dw_z-minority_`y'_`i'*[xb1]f_c9010t_pminorty_dw_z-(${xgamma_short}))
		matrix b = r(b)
		matrix V = r(V)
		scalar povmin2_nox_py`y'_`i'_med_ny = b[1,1]
		scalar povmin2_nox_v`y'_`i' = sqrt(V[1,1])
		scalar povmin2_nox_lci`y'_`i'_med_ny = povmin2_nox_py`y'_`i'_med_ny-invnormal(0.975)*povmin2_nox_v`y'_`i'
		scalar povmin2_nox_uci`y'_`i'_med_ny = povmin2_nox_py`y'_`i'_med_ny+invnormal(0.975)*povmin2_nox_v`y'_`i'
		}
	}

	matrix py`y'_nox_povmin2_med_ny = (scalar(minority_`y'_1),scalar(povmin2_nox_py`y'_1_med_ny),scalar(povmin2_nox_lci`y'_1_med_ny),scalar(povmin2_nox_uci`y'_1_med_ny)\/*
				*/scalar(minority_`y'_2),scalar(povmin2_nox_py`y'_2_med_ny),scalar(povmin2_nox_lci`y'_2_med_ny),scalar(povmin2_nox_uci`y'_2_med_ny)\/*
				*/scalar(minority_`y'_3),scalar(povmin2_nox_py`y'_3_med_ny),scalar(povmin2_nox_lci`y'_3_med_ny),scalar(povmin2_nox_uci`y'_3_med_ny)\/*
				*/scalar(minority_`y'_4),scalar(povmin2_nox_py`y'_4_med_ny),scalar(povmin2_nox_lci`y'_4_med_ny),scalar(povmin2_nox_uci`y'_4_med_ny)\/*
				*/scalar(minority_`y'_5),scalar(povmin2_nox_py`y'_5_med_ny),scalar(povmin2_nox_lci`y'_5_med_ny),scalar(povmin2_nox_uci`y'_5_med_ny)\/*
				*/scalar(minority_`y'_6),scalar(povmin2_nox_py`y'_6_med_ny),scalar(povmin2_nox_lci`y'_6_med_ny),scalar(povmin2_nox_uci`y'_6_med_ny)\/*
				*/scalar(minority_`y'_7),scalar(povmin2_nox_py`y'_7_med_ny),scalar(povmin2_nox_lci`y'_7_med_ny),scalar(povmin2_nox_uci`y'_7_med_ny)\/*
				*/scalar(minority_`y'_8),scalar(povmin2_nox_py`y'_8_med_ny),scalar(povmin2_nox_lci`y'_8_med_ny),scalar(povmin2_nox_uci`y'_8_med_ny)\/*
				*/scalar(minority_`y'_9),scalar(povmin2_nox_py`y'_9_med_ny),scalar(povmin2_nox_lci`y'_9_med_ny),scalar(povmin2_nox_uci`y'_9_med_ny)\/*
				*/scalar(minority_`y'_10),scalar(povmin2_nox_py`y'_10_med_ny),scalar(povmin2_nox_lci`y'_10_med_ny),scalar(povmin2_nox_uci`y'_10_med_ny)\/*
				*/scalar(minority_`y'_11),scalar(povmin2_nox_py`y'_11_med_ny),scalar(povmin2_nox_lci`y'_11_med_ny),scalar(povmin2_nox_uci`y'_11_med_ny))

	matrix colnames py`y'_nox_povmin2_med_ny = povmin2_nox_min`y'_med_ny povmin2_nox_py`y'_med_ny povmin2_nox_lci`y'_med_ny povmin2_nox_uci`y'_med_ny
	svmat py`y'_nox_povmin2_med_ny, names(col)
	
	matrix list py`y'_nox_povmin2_med_ny
}

** ordered probit model 
	
ml model lf myoprobit_lf (xb1:happy_scale012 = $site_covs f_c9010t_perpov_dw_z f_c9010t_pminorty_dw_z, nocons) (cut1:) (cut2:) [pw=$wt]
ml maximize
est store povmin_op_nox

forvalues y = 0/2 {
	if `y'==0{
	margins, at(${x_fix_short} f_c9010t_pminorty_dw_z=.4906762 f_c9010t_perpov_dw_z=(-6(1)4)) expression(normal(predict(equation(cut1))-predict(equation(xb1))))
	}
	else if `y'==1{
	margins, at(${x_fix_short} f_c9010t_pminorty_dw_z=.4906762 f_c9010t_perpov_dw_z=(-6(1)4)) expression((normal(predict(equation(cut2))-predict(equation(xb1))))-(normal(predict(equation(cut1))-predict(equation(xb1)))))
	}
	else if `y'==2{
	margins, at(${x_fix_short} f_c9010t_pminorty_dw_z=.4906762 f_c9010t_perpov_dw_z=(-6(1)4)) expression(1-(normal(predict(equation(cut2))-predict(equation(xb1)))))
	}
	
	marginsplot, ylabel(-0.6(0.3)1.5, grid) title(" ") xtitle(" ") ytitle(" ") xlabel(-6(2)4) recast(line) recastci(rarea) plotopts(lcolor(black) lwidth(medthick)) ciopts(color(black%50)) ///
	addplot(rarea povmin1_nox_lci`y'_med_ny povmin1_nox_uci`y'_med_ny povmin1_nox_pov`y'_med_ny, ylabel(-0.6(0.3)1.5, grid) color(black%20)||line povmin1_nox_py`y'_med_ny ///
	povmin1_nox_pov`y'_med_ny, ylabel(-0.6(0.3)1.5, grid) lcolor(gs5) lpattern(dash) lwidth(medthick)) legend(off)
	graph export ${pathres}prob_y`y'_nox_povmin1_med_ny_${date}.pdf, replace
	}
	
forvalues y = 0/2 {
	if `y'==0{
	margins, at(${x_fix_short} f_c9010t_perpov_dw_z=-.1834148 f_c9010t_pminorty_dw_z=(-6(1)4)) expression(normal(predict(equation(cut1))-predict(equation(xb1))))
	}
	else if `y'==1{
	margins, at(${x_fix_short} f_c9010t_perpov_dw_z=-.1834148 f_c9010t_pminorty_dw_z=(-6(1)4)) expression((normal(predict(equation(cut2))-predict(equation(xb1))))-(normal(predict(equation(cut1))-predict(equation(xb1)))))
	}
	else if `y'==2{
	margins, at(${x_fix_short} f_c9010t_perpov_dw_z=-.1834148 f_c9010t_pminorty_dw_z=(-6(1)4)) expression(1-(normal(predict(equation(cut2))-predict(equation(xb1)))))
	}
	
	marginsplot, ylabel(-0.6(0.3)1.5, grid) title(" ") xtitle(" ") ytitle(" ") xlabel(-6(2)4) recast(line) recastci(rarea) plotopts(lcolor(black) lwidth(medthick)) ciopts(color(black%50)) ///
	addplot(rarea povmin2_nox_lci`y'_med_ny povmin2_nox_uci`y'_med_ny povmin2_nox_min`y'_med_ny, ylabel(-0.6(0.3)1.5, grid) color(black%20)||line povmin2_nox_py`y'_med_ny ///
	povmin2_nox_min`y'_med_ny, ylabel(-0.6(0.3)1.5, grid) lcolor(gs5) lpattern(dash) lwidth(medthick)) legend(off)
	graph export ${pathres}prob_y`y'_nox_povmin2_med_ny_${date}.pdf, replace
	}

** FIGURE 3 

** W = hood poverty, no covariates

** CT model 

ml model lf myoprobit2_lf (xb1:happy_scale012 = $site_covs f_c9010t_perpov_dw_z, nocons) (cut1:) (lndiff:) /*
*/(xb2:f_c9010t_perpov_dw_z = interaction*) (lnsigma:) (rho:) [pw=$wt] , /*
*/diparm(lndiff cut1, f(exp(@1)+@2) d(exp(@1) -exp(@1)+1) label(cut2)) /*
*/diparm(lnsigma, exp label(sigma)) 
ml init xb1:x_f_site_bal=0.1825 x_f_site_bos=0.0528 x_f_site_chi=0.1811 x_f_site_la=0.0674 f_c9010t_perpov_dw_z=-0.1610
ml maximize, difficult
est store pov_t_nox

forvalues y = 0/2 {
	forvalues i = 1/11{
		scalar poverty_`y'_`i' = `i' - 7
		if `y'==0{
		nlcom -[xb1]f_c9010t_perpov_dw_z * normalden([cut1]_cons-poverty_`y'_`i'*[xb1]f_c9010t_perpov_dw_z - (${xgamma_short}) )
		}
		else if `y'==1{
		nlcom [xb1]f_c9010t_perpov_dw_z * (normalden([cut1]_cons-poverty_`y'_`i'*[xb1]f_c9010t_perpov_dw_z - (${xgamma_short}))- normalden(exp([lndiff]_cons)+[cut1]_cons-poverty_`y'_`i'*[xb1]f_c9010t_perpov_dw_z - (${xgamma_short}) ))
		}
		else if `y'==2{
		nlcom [xb1]f_c9010t_perpov_dw_z * normalden(exp([lndiff]_cons)+[cut1]_cons-poverty_`y'_`i'*[xb1]f_c9010t_perpov_dw_z - (${xgamma_short}))
		}
		matrix b = r(b)
		matrix V = r(V)
		scalar me`y'_nox_`i'_med_ny = b[1,1]
		scalar v`y'_nox_`i' = sqrt(V[1,1])
		scalar lci`y'_nox_`i'_med_ny = me`y'_nox_`i'_med_ny-invnormal(0.975)*v`y'_nox_`i'
		scalar uci`y'_nox_`i'_med_ny = me`y'_nox_`i'_med_ny+invnormal(0.975)*v`y'_nox_`i'
	}

	matrix me`y'_nox_pov_med_ny = (scalar(poverty_`y'_1),scalar(me`y'_nox_1_med_ny),scalar(lci`y'_nox_1_med_ny),scalar(uci`y'_nox_1_med_ny)\/*
				*/scalar(poverty_`y'_2),scalar(me`y'_nox_2_med_ny),scalar(lci`y'_nox_2_med_ny),scalar(uci`y'_nox_2_med_ny)\/*
				*/scalar(poverty_`y'_3),scalar(me`y'_nox_3_med_ny),scalar(lci`y'_nox_3_med_ny),scalar(uci`y'_nox_3_med_ny)\/*
				*/scalar(poverty_`y'_4),scalar(me`y'_nox_4_med_ny),scalar(lci`y'_nox_4_med_ny),scalar(uci`y'_nox_4_med_ny)\/*
				*/scalar(poverty_`y'_5),scalar(me`y'_nox_5_med_ny),scalar(lci`y'_nox_5_med_ny),scalar(uci`y'_nox_5_med_ny)\/*
				*/scalar(poverty_`y'_6),scalar(me`y'_nox_6_med_ny),scalar(lci`y'_nox_6_med_ny),scalar(uci`y'_nox_6_med_ny)\/*
				*/scalar(poverty_`y'_7),scalar(me`y'_nox_7_med_ny),scalar(lci`y'_nox_7_med_ny),scalar(uci`y'_nox_7_med_ny)\/*
				*/scalar(poverty_`y'_8),scalar(me`y'_nox_8_med_ny),scalar(lci`y'_nox_8_med_ny),scalar(uci`y'_nox_8_med_ny)\/*
				*/scalar(poverty_`y'_9),scalar(me`y'_nox_9_med_ny),scalar(lci`y'_nox_9_med_ny),scalar(uci`y'_nox_9_med_ny)\/*
				*/scalar(poverty_`y'_10),scalar(me`y'_nox_10_med_ny),scalar(lci`y'_nox_10_med_ny),scalar(uci`y'_nox_10_med_ny)\/*
				*/scalar(poverty_`y'_11),scalar(me`y'_nox_11_med_ny),scalar(lci`y'_nox_11_med_ny),scalar(uci`y'_nox_11_med_ny))

	matrix colnames me`y'_nox_pov_med_ny = mepov`y'_nox_med_ny me`y'_nox_med_ny melci`y'_nox_med_ny meuci`y'_nox_med_ny
	svmat me`y'_nox_pov_med_ny, names(col)
	
	matrix list me`y'_nox_pov_med_ny
}

** ordered probit model 

ml model lf myoprobit_lf (xb1:happy_scale012 = $site_covs f_c9010t_perpov_dw_z, nocons) (cut1:) (cut2:) [pw=$wt]
ml maximize
est store pov_op_nox

forvalues y = 0/2 {
	if `y'==0{
	margins, at(${x_fix_short} f_c9010t_perpov_dw_z=(-6(1)4)) dydx(f_c9010t_perpov_dw_z) expression(normal(predict(equation(cut1))-predict(equation(xb1))))
	}
	else if `y'==1{
	margins, at(${x_fix_short} f_c9010t_perpov_dw_z=(-6(1)4)) dydx(f_c9010t_perpov_dw_z) expression((normal(predict(equation(cut2))-predict(equation(xb1))))-(normal(predict(equation(cut1))-predict(equation(xb1)))))
	}
	else if `y'==2{
	margins, at(${x_fix_short} f_c9010t_perpov_dw_z=(-6(1)4)) dydx(f_c9010t_perpov_dw_z) expression(1-(normal(predict(equation(cut2))-predict(equation(xb1)))))
	}
	
	marginsplot, ylabel(-0.2(0.1)0.2, grid) title(" ") xtitle(" ") ytitle(" ") xlabel(-6(2)4) recast(line) recastci(rarea) plotopts(lcolor(black) lwidth(medthick)) ciopts(color(black%50)) ///
	addplot(rarea melci`y'_nox_med_ny meuci`y'_nox_med_ny mepov`y'_nox_med_ny, ylabel(-0.2(0.1)0.2, grid) color(black%20)||line me`y'_nox_med_ny mepov`y'_nox_med_ny, ///
	ylabel(-0.2(0.1)0.2, grid) lcolor(gs5) lpattern(dash) lwidth(medthick)) legend(off)
	graph export ${pathres}me_y`y'_nox_pov_med_ny_${date}.pdf, replace
	}

** W = hood minority, no covariates

** CT model 

ml model lf myoprobit2_lf (xb1:happy_scale012 = $site_covs f_c9010t_pminorty_dw_z, nocons) (cut1:) (lndiff:) /*
*/(xb2:f_c9010t_pminorty_dw_z = interaction*) (lnsigma:) (rho:) [pw=$wt] , /*
*/diparm(lndiff cut1, f(exp(@1)+@2) d(exp(@1) -exp(@1)+1) label(cut2)) /*
*/diparm(lnsigma, exp label(sigma)) 
ml maximize
est store min_t_nox

forvalues y = 0/2{
	forvalues i = 1/11{
		scalar minority_`y'_`i' = `i' - 7
		if `y'==0{
		nlcom -[xb1]f_c9010t_pminorty_dw_z * normalden([cut1]_cons-minority_`y'_`i'*[xb1]f_c9010t_pminorty_dw_z - (${xgamma_short}) )
		}
		else if `y'==1{
		nlcom [xb1]f_c9010t_pminorty_dw_z * (normalden([cut1]_cons-minority_`y'_`i'*[xb1]f_c9010t_pminorty_dw_z - (${xgamma_short}) )- /*
				*/normalden(exp([lndiff]_cons)+[cut1]_cons-minority_`y'_`i'*[xb1]f_c9010t_pminorty_dw_z - (${xgamma_short}) ))
		}
		else if `y'==2{
		nlcom [xb1]f_c9010t_pminorty_dw_z * normalden(exp([lndiff]_cons)+[cut1]_cons-minority_`y'_`i'*[xb1]f_c9010t_pminorty_dw_z - (${xgamma_short}) )
		}
		matrix b = r(b)
		matrix V = r(V)
		scalar me`y'_nox_`i'_med_ny = b[1,1]
		scalar v`y'_nox_`i' = sqrt(V[1,1])
		scalar lci`y'_nox_`i'_med_ny = me`y'_nox_`i'_med_ny-invnormal(0.975)*v`y'_nox_`i'
		scalar uci`y'_nox_`i'_med_ny = me`y'_nox_`i'_med_ny+invnormal(0.975)*v`y'_nox_`i'
	}

	matrix me`y'_nox_min_med_ny = (scalar(minority_`y'_1),scalar(me`y'_nox_1_med_ny),scalar(lci`y'_nox_1_med_ny),scalar(uci`y'_nox_1_med_ny)\/*
				*/scalar(minority_`y'_2),scalar(me`y'_nox_2_med_ny),scalar(lci`y'_nox_2_med_ny),scalar(uci`y'_nox_2_med_ny)\/*
				*/scalar(minority_`y'_3),scalar(me`y'_nox_3_med_ny),scalar(lci`y'_nox_3_med_ny),scalar(uci`y'_nox_3_med_ny)\/*
				*/scalar(minority_`y'_4),scalar(me`y'_nox_4_med_ny),scalar(lci`y'_nox_4_med_ny),scalar(uci`y'_nox_4_med_ny)\/*
				*/scalar(minority_`y'_5),scalar(me`y'_nox_5_med_ny),scalar(lci`y'_nox_5_med_ny),scalar(uci`y'_nox_5_med_ny)\/*
				*/scalar(minority_`y'_6),scalar(me`y'_nox_6_med_ny),scalar(lci`y'_nox_6_med_ny),scalar(uci`y'_nox_6_med_ny)\/*
				*/scalar(minority_`y'_7),scalar(me`y'_nox_7_med_ny),scalar(lci`y'_nox_7_med_ny),scalar(uci`y'_nox_7_med_ny)\/*
				*/scalar(minority_`y'_8),scalar(me`y'_nox_8_med_ny),scalar(lci`y'_nox_8_med_ny),scalar(uci`y'_nox_8_med_ny)\/*
				*/scalar(minority_`y'_9),scalar(me`y'_nox_9_med_ny),scalar(lci`y'_nox_9_med_ny),scalar(uci`y'_nox_9_med_ny)\/*
				*/scalar(minority_`y'_10),scalar(me`y'_nox_10_med_ny),scalar(lci`y'_nox_10_med_ny),scalar(uci`y'_nox_10_med_ny)\/*
				*/scalar(minority_`y'_11),scalar(me`y'_nox_11_med_ny),scalar(lci`y'_nox_11_med_ny),scalar(uci`y'_nox_11_med_ny))

	matrix colnames me`y'_nox_min_med_ny = memin`y'_nox_med_ny minme`y'_nox_med_ny minmelci`y'_nox_med_ny minmeuci`y'_nox_med_ny
	svmat me`y'_nox_min_med_ny, names(col)
	
	matrix list me`y'_nox_min_med_ny
}

** ordered probit model 

ml model lf myoprobit_lf (xb1:happy_scale012 = $site_covs f_c9010t_pminorty_dw_z, nocons) (cut1:) (cut2:) [pw=$wt]
ml maximize
est store min_op_nox
	
forvalues y = 0/2 {
	if `y'==0{
	margins, at(${x_fix_short} f_c9010t_pminorty_dw_z=(-6(1)4)) dydx(f_c9010t_pminorty_dw_z) expression(normal(predict(equation(cut1))-predict(equation(xb1))))
	}
	else if `y'==1{
	margins, at(${x_fix_short} f_c9010t_pminorty_dw_z=(-6(1)4)) dydx(f_c9010t_pminorty_dw_z) expression((normal(predict(equation(cut2))-predict(equation(xb1))))-(normal(predict(equation(cut1))-predict(equation(xb1)))))
	}
	else if `y'==2{
	margins, at(${x_fix_short} f_c9010t_pminorty_dw_z=(-6(1)4)) dydx(f_c9010t_pminorty_dw_z) expression(1-(normal(predict(equation(cut2))-predict(equation(xb1)))))
	}
	
	marginsplot, ylabel(-0.2(0.1)0.2, grid) title(" ") xtitle(" ") ytitle(" ") xlabel(-6(2)4) recast(line) recastci(rarea) plotopts(lcolor(black) lwidth(medthick)) ciopts(color(black%50)) ///
	addplot(rarea minmelci`y'_nox_med_ny minmeuci`y'_nox_med_ny memin`y'_nox_med_ny, ylabel(-0.2(0.1)0.2, grid) color(black%20)||line minme`y'_nox_med_ny memin`y'_nox_med_ny, ///
	ylabel(-0.2(0.1)0.2, grid) lcolor(gs5) lpattern(dash) lwidth(medthick)) legend(off)
	graph export ${pathres}me_y`y'_nox_min_med_ny_${date}.pdf, replace
	}

** FIGURE 4 

** W = (hood poverty, hood minority), no covariates

** CT model 

ml model lf myoprobit3_lf (xb1:happy_scale012 = $site_covs f_c9010t_perpov_dw_z f_c9010t_pminorty_dw_z, nocons) (cut1:) (lndiff:) /*
						*/(xb2:f_c9010t_perpov_dw_z = interaction*) (lnvar_v1:) (rho_uv1:) /*
						*/(xb3:f_c9010t_pminorty_dw_z = interaction*) (lnvar_v2:) (rho_uv2:) (cov_v1v2:) [pw=$wt], /*
						*/diparm(lndiff cut1, f(exp(@1)+@2) d(exp(@1) -exp(@1)+1) label(cut2)) /*
						*/diparm(lnvar_v1, exp label(var_v1)) /*
						*/diparm(lnvar_v2, exp label(var_v2)) 
ml init xb1:x_f_site_bal=0.3039 x_f_site_bos=0.3386 x_f_site_chi=0.1527 x_f_site_la=0.0639 f_c9010t_perpov_dw_z=-0.2998 f_c9010t_pminorty_dw_z=0.3151 cut1:_cons=-0.8165
ml maximize, difficult
est store povmin_t_nox

forvalues y = 0/2{
	forvalues i = 1/11{
		scalar poverty_`y'_`i' = `i' - 7
		if `y'==0{
		nlcom -[xb1]f_c9010t_perpov_dw_z * normalden([cut1]_cons-poverty_`y'_`i'*[xb1]f_c9010t_perpov_dw_z-f_c9010t_pminorty_dw_z_med_nox*[xb1]f_c9010t_pminorty_dw_z-(${xgamma_short}))
		matrix b = r(b)
		matrix V = r(V)
		scalar povmin1_nox_me`y'_`i'_med_ny = b[1,1]
		scalar povmin1_nox_v`y'_`i' = sqrt(V[1,1])
		scalar povmin1_nox_lci`y'_`i'_med_ny = povmin1_nox_me`y'_`i'_med_ny-invnormal(0.975)*povmin1_nox_v`y'_`i'
		scalar povmin1_nox_uci`y'_`i'_med_ny = povmin1_nox_me`y'_`i'_med_ny+invnormal(0.975)*povmin1_nox_v`y'_`i'
		}
		else if `y'==1{
		nlcom [xb1]f_c9010t_perpov_dw_z * normalden([cut1]_cons-poverty_`y'_`i'*[xb1]f_c9010t_perpov_dw_z-f_c9010t_pminorty_dw_z_med_nox*[xb1]f_c9010t_pminorty_dw_z-(${xgamma_short}))- /*
				*/[xb1]f_c9010t_perpov_dw_z * normalden(exp([lndiff]_cons)+[cut1]_cons-poverty_`y'_`i'*[xb1]f_c9010t_perpov_dw_z-f_c9010t_pminorty_dw_z_med_nox*[xb1]f_c9010t_pminorty_dw_z-(${xgamma_short}))
		matrix b = r(b)
		matrix V = r(V)
		scalar povmin1_nox_me`y'_`i'_med_ny = b[1,1]
		scalar povmin1_nox_v`y'_`i' = sqrt(V[1,1])
		scalar povmin1_nox_lci`y'_`i'_med_ny = povmin1_nox_me`y'_`i'_med_ny-invnormal(0.975)*povmin1_nox_v`y'_`i'
		scalar povmin1_nox_uci`y'_`i'_med_ny = povmin1_nox_me`y'_`i'_med_ny+invnormal(0.975)*povmin1_nox_v`y'_`i'
		}
		else if `y'==2{
		nlcom [xb1]f_c9010t_perpov_dw_z * normalden(exp([lndiff]_cons)+[cut1]_cons-poverty_`y'_`i'*[xb1]f_c9010t_perpov_dw_z-f_c9010t_pminorty_dw_z_med_nox*[xb1]f_c9010t_pminorty_dw_z-(${xgamma_short}))
		matrix b = r(b)
		matrix V = r(V)
		scalar povmin1_nox_me`y'_`i'_med_ny = b[1,1]
		scalar povmin1_nox_v`y'_`i' = sqrt(V[1,1])
		scalar povmin1_nox_lci`y'_`i'_med_ny = povmin1_nox_me`y'_`i'_med_ny-invnormal(0.975)*povmin1_nox_v`y'_`i'
		scalar povmin1_nox_uci`y'_`i'_med_ny = povmin1_nox_me`y'_`i'_med_ny+invnormal(0.975)*povmin1_nox_v`y'_`i'
		}
	}

	matrix me`y'_povmin1_nox_med_ny = (scalar(poverty_`y'_1),scalar(povmin1_nox_me`y'_1_med_ny),scalar(povmin1_nox_lci`y'_1_med_ny),scalar(povmin1_nox_uci`y'_1_med_ny)\/*
				*/scalar(poverty_`y'_2),scalar(povmin1_nox_me`y'_2_med_ny),scalar(povmin1_nox_lci`y'_2_med_ny),scalar(povmin1_nox_uci`y'_2_med_ny)\/*
				*/scalar(poverty_`y'_3),scalar(povmin1_nox_me`y'_3_med_ny),scalar(povmin1_nox_lci`y'_3_med_ny),scalar(povmin1_nox_uci`y'_3_med_ny)\/*
				*/scalar(poverty_`y'_4),scalar(povmin1_nox_me`y'_4_med_ny),scalar(povmin1_nox_lci`y'_4_med_ny),scalar(povmin1_nox_uci`y'_4_med_ny)\/*
				*/scalar(poverty_`y'_5),scalar(povmin1_nox_me`y'_5_med_ny),scalar(povmin1_nox_lci`y'_5_med_ny),scalar(povmin1_nox_uci`y'_5_med_ny)\/*
				*/scalar(poverty_`y'_6),scalar(povmin1_nox_me`y'_6_med_ny),scalar(povmin1_nox_lci`y'_6_med_ny),scalar(povmin1_nox_uci`y'_6_med_ny)\/*
				*/scalar(poverty_`y'_7),scalar(povmin1_nox_me`y'_7_med_ny),scalar(povmin1_nox_lci`y'_7_med_ny),scalar(povmin1_nox_uci`y'_7_med_ny)\/*
				*/scalar(poverty_`y'_8),scalar(povmin1_nox_me`y'_8_med_ny),scalar(povmin1_nox_lci`y'_8_med_ny),scalar(povmin1_nox_uci`y'_8_med_ny)\/*
				*/scalar(poverty_`y'_9),scalar(povmin1_nox_me`y'_9_med_ny),scalar(povmin1_nox_lci`y'_9_med_ny),scalar(povmin1_nox_uci`y'_9_med_ny)\/*
				*/scalar(poverty_`y'_10),scalar(povmin1_nox_me`y'_10_med_ny),scalar(povmin1_nox_lci`y'_10_med_ny),scalar(povmin1_nox_uci`y'_10_med_ny)\/*
				*/scalar(poverty_`y'_11),scalar(povmin1_nox_me`y'_11_med_ny),scalar(povmin1_nox_lci`y'_11_med_ny),scalar(povmin1_nox_uci`y'_11_med_ny))

	matrix colnames me`y'_povmin1_nox_med_ny = povmin1_nox_mepov`y'_med_ny povmin1_nox_me`y'_med_ny povmin1_nox_melci`y'_med_ny povmin1_nox_meuci`y'_med_ny
	svmat me`y'_povmin1_nox_med_ny, names(col)

	matrix list me`y'_povmin1_nox_med_ny
}

forvalues y = 0/2{
	forvalues i = 1/11{
		scalar minority_`y'_`i' = `i' - 7
		if `y'==0{
		nlcom -[xb1]f_c9010t_pminorty_dw_z * normalden([cut1]_cons-f_c9010t_perpov_dw_z_med_nox*[xb1]f_c9010t_perpov_dw_z-minority_`y'_`i'*[xb1]f_c9010t_pminorty_dw_z-(${xgamma_short}))
		matrix b = r(b)
		matrix V = r(V)
		scalar povmin2_nox_me`y'_`i'_med_ny = b[1,1]
		scalar povmin2_nox_v`y'_`i' = sqrt(V[1,1])
		scalar povmin2_nox_lci`y'_`i'_med_ny = povmin2_nox_me`y'_`i'_med_ny-invnormal(0.975)*povmin2_nox_v`y'_`i'
		scalar povmin2_nox_uci`y'_`i'_med_ny = povmin2_nox_me`y'_`i'_med_ny+invnormal(0.975)*povmin2_nox_v`y'_`i'
		}
		else if `y'==1{
		nlcom [xb1]f_c9010t_pminorty_dw_z * normalden([cut1]_cons-f_c9010t_perpov_dw_z_med_nox*[xb1]f_c9010t_perpov_dw_z-minority_`y'_`i'*[xb1]f_c9010t_pminorty_dw_z-(${xgamma_short})) - /*
				*/[xb1]f_c9010t_pminorty_dw_z * normalden(exp([lndiff]_cons)+[cut1]_cons-f_c9010t_perpov_dw_z_med_nox*[xb1]f_c9010t_perpov_dw_z-minority_`y'_`i'*[xb1]f_c9010t_pminorty_dw_z-(${xgamma_short}))
		matrix b = r(b)
		matrix V = r(V)
		scalar povmin2_nox_me`y'_`i'_med_ny = b[1,1]
		scalar povmin2_nox_v`y'_`i' = sqrt(V[1,1])
		scalar povmin2_nox_lci`y'_`i'_med_ny = povmin2_nox_me`y'_`i'_med_ny-invnormal(0.975)*povmin2_nox_v`y'_`i'
		scalar povmin2_nox_uci`y'_`i'_med_ny = povmin2_nox_me`y'_`i'_med_ny+invnormal(0.975)*povmin2_nox_v`y'_`i'
		}
		else if `y'==2{
		nlcom [xb1]f_c9010t_pminorty_dw_z * normalden(exp([lndiff]_cons)+[cut1]_cons-f_c9010t_perpov_dw_z_med_nox*[xb1]f_c9010t_perpov_dw_z-minority_`y'_`i'*[xb1]f_c9010t_pminorty_dw_z-(${xgamma_short}))
		matrix b = r(b)
		matrix V = r(V)
		scalar povmin2_nox_me`y'_`i'_med_ny = b[1,1]
		scalar povmin2_nox_v`y'_`i' = sqrt(V[1,1])
		scalar povmin2_nox_lci`y'_`i'_med_ny = povmin2_nox_me`y'_`i'_med_ny-invnormal(0.975)*povmin2_nox_v`y'_`i'
		scalar povmin2_nox_uci`y'_`i'_med_ny = povmin2_nox_me`y'_`i'_med_ny+invnormal(0.975)*povmin2_nox_v`y'_`i'
		}
	}

	matrix me`y'_povmin2_nox_med_ny = (scalar(minority_`y'_1),scalar(povmin2_nox_me`y'_1_med_ny),scalar(povmin2_nox_lci`y'_1_med_ny),scalar(povmin2_nox_uci`y'_1_med_ny)\/*
				*/scalar(minority_`y'_2),scalar(povmin2_nox_me`y'_2_med_ny),scalar(povmin2_nox_lci`y'_2_med_ny),scalar(povmin2_nox_uci`y'_2_med_ny)\/*
				*/scalar(minority_`y'_3),scalar(povmin2_nox_me`y'_3_med_ny),scalar(povmin2_nox_lci`y'_3_med_ny),scalar(povmin2_nox_uci`y'_3_med_ny)\/*
				*/scalar(minority_`y'_4),scalar(povmin2_nox_me`y'_4_med_ny),scalar(povmin2_nox_lci`y'_4_med_ny),scalar(povmin2_nox_uci`y'_4_med_ny)\/*
				*/scalar(minority_`y'_5),scalar(povmin2_nox_me`y'_5_med_ny),scalar(povmin2_nox_lci`y'_5_med_ny),scalar(povmin2_nox_uci`y'_5_med_ny)\/*
				*/scalar(minority_`y'_6),scalar(povmin2_nox_me`y'_6_med_ny),scalar(povmin2_nox_lci`y'_6_med_ny),scalar(povmin2_nox_uci`y'_6_med_ny)\/*
				*/scalar(minority_`y'_7),scalar(povmin2_nox_me`y'_7_med_ny),scalar(povmin2_nox_lci`y'_7_med_ny),scalar(povmin2_nox_uci`y'_7_med_ny)\/*
				*/scalar(minority_`y'_8),scalar(povmin2_nox_me`y'_8_med_ny),scalar(povmin2_nox_lci`y'_8_med_ny),scalar(povmin2_nox_uci`y'_8_med_ny)\/*
				*/scalar(minority_`y'_9),scalar(povmin2_nox_me`y'_9_med_ny),scalar(povmin2_nox_lci`y'_9_med_ny),scalar(povmin2_nox_uci`y'_9_med_ny)\/*
				*/scalar(minority_`y'_10),scalar(povmin2_nox_me`y'_10_med_ny),scalar(povmin2_nox_lci`y'_10_med_ny),scalar(povmin2_nox_uci`y'_10_med_ny)\/*
				*/scalar(minority_`y'_11),scalar(povmin2_nox_me`y'_11_med_ny),scalar(povmin2_nox_lci`y'_11_med_ny),scalar(povmin2_nox_uci`y'_11_med_ny))

	matrix colnames me`y'_povmin2_nox_med_ny = povmin2_nox_memin`y'_med_ny povmin2_nox_me`y'_med_ny povmin2_nox_melci`y'_med_ny povmin2_nox_meuci`y'_med_ny
	svmat me`y'_povmin2_nox_med_ny, names(col)

	matrix list me`y'_povmin2_nox_med_ny
}

** ordered probit 
	
ml model lf myoprobit_lf (xb1:happy_scale012 = $site_covs f_c9010t_perpov_dw_z f_c9010t_pminorty_dw_z, nocons) (cut1:) (cut2:) [pw=$wt]
ml maximize
est store povmin_op_nox

forvalues y = 0/2 {
	if `y'==0{
	margins, at(${x_fix_short} f_c9010t_pminorty_dw_z=.4906762 f_c9010t_perpov_dw_z=(-6(1)4)) dydx(f_c9010t_perpov_dw_z) expression(normal(predict(equation(cut1))-predict(equation(xb1))))
	}
	else if `y'==1{
	margins, at(${x_fix_short} f_c9010t_pminorty_dw_z=.4906762 f_c9010t_perpov_dw_z=(-6(1)4)) dydx(f_c9010t_perpov_dw_z) expression((normal(predict(equation(cut2))-predict(equation(xb1))))-(normal(predict(equation(cut1))-predict(equation(xb1)))))
	}
	else if `y'==2{
	margins, at(${x_fix_short} f_c9010t_pminorty_dw_z=.4906762 f_c9010t_perpov_dw_z=(-6(1)4)) dydx(f_c9010t_perpov_dw_z) expression(1-(normal(predict(equation(cut2))-predict(equation(xb1)))))
	}
	
	marginsplot, ylabel(-0.4(0.2)0.4, grid) title(" ") xtitle(" ") ytitle(" ") xlabel(-6(2)4) recast(line) recastci(rarea) plotopts(lcolor(black) lwidth(medthick)) ciopts(color(black%50)) ///
	addplot(rarea povmin1_nox_melci`y'_med_ny povmin1_nox_meuci`y'_med_ny povmin1_nox_mepov`y'_med_ny, ylabel(-0.4(0.2)0.4, grid) color(black%20)||line povmin1_nox_me`y'_med_ny ///
	povmin1_nox_mepov`y'_med_ny, ylabel(-0.4(0.2)0.4, grid) lcolor(gs5) lpattern(dash) lwidth(medthick)) legend(off)
	graph export ${pathres}me_y`y'_nox_povmin1_med_ny_${date}.pdf, replace
	}
	
forvalues y = 0/2 {
	if `y'==0{
	margins, at(${x_fix_short} f_c9010t_perpov_dw_z=-.1834148 f_c9010t_pminorty_dw_z=(-6(1)4)) dydx(f_c9010t_pminorty_dw_z) expression(normal(predict(equation(cut1))-predict(equation(xb1))))
	}
	else if `y'==1{
	margins, at(${x_fix_short} f_c9010t_perpov_dw_z=-.1834148 f_c9010t_pminorty_dw_z=(-6(1)4)) dydx(f_c9010t_pminorty_dw_z) expression((normal(predict(equation(cut2))-predict(equation(xb1))))-(normal(predict(equation(cut1))-predict(equation(xb1)))))
	}
	else if `y'==2{
	margins, at(${x_fix_short} f_c9010t_perpov_dw_z=-.1834148 f_c9010t_pminorty_dw_z=(-6(1)4)) dydx(f_c9010t_pminorty_dw_z) expression(1-(normal(predict(equation(cut2))-predict(equation(xb1)))))
	}
	
	marginsplot, ylabel(-0.4(0.2)0.4, grid) title(" ") xtitle(" ") ytitle(" ") xlabel(-6(2)4) recast(line) recastci(rarea) plotopts(lcolor(black) lwidth(medthick)) ciopts(color(black%50)) ///
	addplot(rarea povmin2_nox_melci`y'_med_ny povmin2_nox_meuci`y'_med_ny povmin2_nox_memin`y'_med_ny, ylabel(-0.4(0.2)0.4, grid) color(black%20)||line povmin2_nox_me`y'_med_ny ///
	povmin2_nox_memin`y'_med_ny, ylabel(-0.4(0.2)0.4, grid) lcolor(gs5) lpattern(dash) lwidth(medthick)) legend(off)
	graph export ${pathres}me_y`y'_nox_povmin2_med_ny_${date}.pdf, replace
	}

scalar drop _all
matrix drop _all
estimates drop _all
